	use "${clean_data}/panel_physicians_clean.dta", clear
		
	encode(cms_spec_code), gen(cms_spec_num)
	encode(stfips), gen(stfips_num)
	
	cap log close
	log using "${output}/descriptives/time_series_regs.txt", replace text 	
	
	replace ptotinc = ptotinc/1000

	reg ptotinc ib2017.year i.age female i.cms_spec_num i.stfips_num, robust 
	parmest, saving("${output}/descriptives/time_series_adjusted_ptotinc.dta", replace) es(N)
	gen adjusted_sample=(e(sample)==1)

	reg top1 ib2017.year i.age female i.cms_spec_num i.stfips_num if adjusted_sample==1, robust 
	parmest, saving("${output}/descriptives/time_series_adjusted_top1.dta", replace) es(N)

	reg ptotinc ib2017.year if adjusted_sample==1
	parmest, saving("${output}/descriptives/time_series_notadjusted_ptotinc.dta", replace) es(N)
	
	reg top1 ib2017.year if adjusted_sample==1
	parmest, saving("${output}/descriptives/time_series_notadjusted_top1.dta", replace)	es(N)
	
	reg top1 ib2017.year if einsize==.
	parmest, saving("${output}/descriptives/time_series_notadjusted_einmissing.dta", replace)	es(N)
	
	log close
	
	* Plot time series

	use "${output}/descriptives/time_series_notadjusted_ptotinc.dta", clear	
	
	keep parm estimate es_1
	rename (estimate es_1) (estimate_notadj_ptotinc N_notadj_ptotinc_UR)
	merge 1:1 parm using  "${output}/descriptives/time_series_adjusted_ptotinc.dta", keep(master match) keepusing(estimate es_1) nogen
	rename (estimate es_1) (estimate_adj_ptotinc N_adj_ptotinc_UR)
	
	merge 1:1 parm using  "${output}/descriptives/time_series_notadjusted_top1.dta", keep(master match) keepusing(estimate es_1) nogen
	rename (estimate es_1) (estimate_notadj_top1 N_notadj_top1_UR)
	
	merge 1:1 parm using  "${output}/descriptives/time_series_adjusted_top1.dta", keep(master match) keepusing(estimate es_1) nogen
	rename (estimate es_1) (estimate_adj_top1 N_adj_top1_UR)

	gen mean2017_ptotoinc = estimate_notadj_ptotinc[14]
	gen mean2017_top1 = estimate_notadj_top1[14]
	
	drop if parm == "_cons"
	
	gen notadj_ptotinc=estimate_notadj_ptotinc+mean2017_ptotoinc
	gen adj_ptotinc=estimate_adj_ptotinc+mean2017_ptotoinc
	gen notadj_top1=estimate_notadj_top1+mean2017_top1
	gen adj_top1=estimate_adj_top1+mean2017_top1
		
	gen year=substr(parm, 1,4)
	destring year, replace
	
	* Export dataset
	
	keep notadj_ptotinc adj_ptotinc notadj_top1 adj_top1 year N*
	drbest_docinc notadj_ptotinc adj_ptotinc notadj_top1 adj_top1
	order N* , last
	export delimited using "$mypath/intermediate_csv/desc07-time_series_earnings.csv", replace dataf delim(tab)  
